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Large Eddy Simulations and Direct Numerical 
Simulations of High Speed Turbulent Reacting 

Flows 


P. Givi, C.K. Madnia, C .J. Steinberger, S.H. Frankel and T.J. Vidoni 
Department of Mechanical and Aerospace Engineering 
State University of New York at Buffalo 
Buffalo, New York 14260 


Abstract 

This main objective of this research is to extend the boundaries within which 
Large Eddy Simulations (LES) and Direct Numerical Simulations (DNS) can be applied 
in computational analyses of high speed reacting flows. In the efforts related to LES, 
we have been concerned with developing reliable subgrid closures for modeling of 
the fluctuation correlations of scalar quantities in reacting turbulent flows. In the 
work on DNS, we have focused our attention to further investigate the effects of 
exothermicity in compressible turbulent flows. In our previous work, in the first year 
of this research, we have considered only "simple" flows. Currently, we are in. the 
process of extending our analyses for the purpose of modeling more practical flows 
of current interest at NASA Langley Research Center. 

This report provides a summary of our accomplishments during the third six 
months of the research supported under Grant NAG 1-1122 sponsored by NASA 
Langley Research Center administrated by Dr. J. Phil Drummond of the Theoretical 
Flow Physics Branch. 
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1 Introduction 


During the second year of activities under this research, we have been primarily involved 
with the continuations of our tasks in the same two directions as those followed in the 
first year. Namely: (1) development of subgrid closures for LES of compressible reacting 
flows, and (2) understanding the mechanisms of mixing and reaction in high speed 
combustion. In the efforts related to the first task, we have undertaken the simulations of 
homogeneous reacting turbulence, whereas in the second task, DNS of a reacting planar 
jet is the subject of main focus. Below, a summary of our efforts during the first 6 months 
of this research in the second year is provided and detail descriptions are furnished in 
Appendix I and Appendix II. 


1.1 Large Eddy Simulations of Compressible Reacting Flows 

Our major goal in this effort is to initiate a program to extend the present capabilities 
of LES for the treatment of chemically reacting flows. In the efforts to date, we have 
been primarily concerned with a priori analysis of subgrid fluctuation in a compressible 
homogeneous flows. This analysis is mainly involved with constructing the shapes of 
the simulated Probability Density Functions (PDF's) within the subgrid. Based on our 
efforts thus far, we have been able to demonstrate that in the context of single-point 
PDF formulation, the mapping closure of Kraichnan [1] provides the most satisfactory 
results. The main deficiency of this mapping closure (or any other single-point PDF 
model) is associated with its incapability in predicting the frequency (or length) scale of 
PDF evolution. At this stage, unfortunately, the extension of the mapping closure for PDF 
formulation at two-point is very difficult. Among the currently available alternatives, the 
EDQNM spectral closure [2] provides a reasonable remedy to overcom this deficiency, 
but at the present stage cannot be yet implemented for the modeling of compressible 
turbulent flows. 

With the utilization of the mapping closure, the limiting rate of fuel consumption in 
homogeneous flows can be predicted by a closed form mathematical relation. The results 
obtained by this relation compares very well with DNS data. In the firs year, our 
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results were valid only for stoichiometric mixtures. In the second year, we have extended 
our analysis for the modeling of non-unity equivalence ratio mixtures. To the best of 
our knowledge, this relation is more satisfactory than any other suggested correlations 
(empirical or theoretical) available in the literature within the past thirty years! 

In addition to these conclusions, we have also been able to demonstrate that the ^-density 
[3] does reasonably well in approximating the single point PDF distribution of a conserved 
scalar property, if the evolution of its variance is known a priori. With the implementation 
of this density, it is also possible to develop a closed form algebraic solution for the 
limiting rate of reactant conversion under all equivalence ratios. However, despite its 
simplicity the ^-density is only applicable for description of "simple" distributions. It is 
not applicable for trimodal distributions, and it can not predict the subsequent evolution 
of the PDF, if the distribution is not of /?-type initially. Furthermore, it can not be extended 
for predictions of non-equilibrium reactions and there are no firm mathematical proofs to 
justify its use. 

At this stage, we have not yet applied the mapping closure for the modeling of practical 
flows. However, the model based on the /0-density has proven somewhat useful for the 
purpose of turbulence modeling in free shear flows. Appendix I provides a summary of 
our activities in this regard. The extension of our work for the purpose of subgrid closures 
in such flows, is the subject of our current investigation. 


1.2 Direct Numerical Simulations of an Exothermic Reacting Planar Jet 

Our efforts in this work have been focused on elucidating the role of exothermicity in 
an unpremixed reacting planar jet flame. This work has been mainly motivated due to 
the interest of NASA in understanding the global and detail mechanisms of mixing and 
chemical reactions in high speed jet flames. 

In this work, results are obtained of DNS of a compressible, spatially developing, 
forced planar jet under the influence of a finite rate chemical reaction of the type 
F + O — > Product + Heat. Our results suggest that in the setting of a "turbulent" flame, 
the effect of the heat liberated by the chemical reaction is to increase the rate of reactant 
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conversion. This finding is different from the conclusions of earlier DNS results and 
laboratory investigations which indicate a suppressed chemical reaction with increasing 
heat release. However, this conclusion is established only in the context of the simplified 
kinetics mechanism adopted here, and its extent of applicability for practical flames 
requires laboratory measurements. 

We are presently at a preliminary stage of our work on this issue. A summary of our 
findings to date is provided in Appendix II. 


2 Publications and Honors 

The following publications have resulted from our efforts during the first phase of the 
second year of this research: 

1. S.H. Frankel, C.K. Madnia, and P. Givi, "Modeling of the Reactant Conversion Rate in 
a Turbulent Shear Flow," Chem. Eng. Comm., in press (1992). 

2. C.J. Steinberger, "Model Free Simulations of a High Speed Reacting Mixing Layer", 
AIAA Paper 92-0257 (1992). 

Mr. Craig Steinberger was invited to participate at the AIAA 30th Aerospace Sciences 
Meeting, held in Reno, Nevada, to present his AIAA paper listed above. 
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3 Appendix I 


This appendix has been prepared by Mr. Steven Frankel and Dr. Cyrus Madnia. Dr. 
Madnia is currently being partially supported by this grant. 


Problem Description 


The problem under consideration is that of a spatially evolving mixing layer containing 
initially segregated reactants A and B in the two free streams. Species A is introduced 
into the upper, high-speed stream and species B enters on the lower, low-speed side. The 
chemical reaction between the two species is assumed to be single step, irreversible and 
infinitely fast. The magnitude of the Mach number is assumed negligible. Therefore, the 
flow is considered incompressible. The two species are assumed thermodynamically sim- 
ilar with identical diffusion coefficients. Within the framework of these approximations, 
the flow is mathematically described by the conservation equations of mass, momentum, 
and a species conservation equation for the trace of a conserved Shvab-Zeldovich variable 
J, which characterizes the compositional structure of the flow. These equations are 
parameterized by the Reynolds number, the Peclet number, and the velocity ratio across 
the layer. 

The computational package employed in the DNS is based on a hybrid pseudospectral- 
spectral element algorithm developed in [4, 5]. The pseudospectral routine, which is 
based on Fourier collocation, is used in the cross stream direction together with free slip 
boundary conditions. The spectral-element discretization is employed in the streamwise 
direction. This involves the decomposition of the domain in this direction into an 
ensemble of macro finite elements. Within each of these elements, the dependent flow 
variables are approximated spectrally by means of Tchebysheff polynomials of the first 
kind. The hybrid procedure employed in this way is very attractive in that it combines 
the accuracy of spectral discretization with the versatility of finite element methods. 
Therefore, it is convenient for accurate simulations of the complex spatially developing 
flow under consideration here. For a detailed description of this hybrid method for DNS 
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of turbulent shear flows, the reader is referred to [5]. 


The flow field is initialized at the inflow by a hyperbolic tangent mean velocity distribution, 
upon which low amplitude perturbations are superimposed in order to promote the 
formation of coherent vortices. The initialization of the Shvab-Zeldovich parameter at the 
inflow is such that it takes on the values of zero and one in the streams containing B and 
A, respectively. The disturbances on the mean flow correspond to the most unstable mode 
of the hyperbolic tangent profile [6] and its first two subharmonics. The amplitude of the 
fluctuating velocity is set equal to 6% of the mean velocity. This amplitude corresponds 
to that measured in typical shear layer experiments. The magnitudes of the phase shifts 
between the modes of the instability waves are randomly selected from a random seed 
with a top hat PDF of zero mean and a specified variance. The implementation of 
these phase shifts is the only mechanism to introduce randomness into an otherwise 
deterministic simulation. This is to partially mimic the random "commotions" which 
are present in the "universe," into an isolated two-dimensional simulation. At the 
outflow, a weak condition of zero second derivative is applied for all the dependent 
variables. This boundary condition cannot be justified in a rigorous mathematical sense 
and can be specified only in an ad hoc manner. This was done in order to facilitate 
implementation of the Dirichlet boundary conditions in the finite element procedure 
employed here. Moreover, the results of extensive numerical tests showed that the effects 
of the approximate boundary condition are confined to within the last computational 
element. This is consistent with that found previously in [7]. Therefore, the solution in 
the last computational domain can be ignored. 

With the solution of the unsteady transport equations, the DNS data are extremely 
useful in visualizing the instantaneous behavior of the flow. The results of these direct 
simulations can also provide an assessment of the statistical behavior of the flow. This is 
due to the availability of flow information at all the computational grid points and at all 
times. Therefore, statistical sampling can be done quite easily by ensemble averaging of 
the instantaneous data gathered over many realizations. 

In our simulations the instantaneous data of the Shvab-Zeldovich variable, obtained from 
the DNS, are used for statistical analysis. A total of 3200 realizations were employed in 
the sampling. With this ensemble, the magnitudes of the mean and the variance of the 
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Shvab-Zeldovich variable are calculated at all the grid points. These statistical quantities 
are of primary interest for turbulence modeling. Based on the knowledge of these first 
two moments, a Beta density is assumed for the PDF of the Shvab-Zeldovich variable. The 
ensemble mean values of the product concentration and the unmixedness are calculated 
subsequently via this density. These modeled quantities are then compared with those 
obtained directly from the DNS data. Due to the spatial inhomogeneity of the flow the 
mixture is not stoichiometric at all the grid points. Therefore, the Beta density is usually 
asymmetric at majority of these points. 


Presentation of Results 

Computations were performed in a domain of size 135 x 345, where 5 denotes the 
vorticity thickness at the inlet of the flow. There are 65 Fourier modes in the cross-stream 
direction and 42 finite elements are used for the streamwise discretization. A fifth order 
Tchebysheff polynomial is employed to approximate the variables within each element. 
This discretization is equivalent to, atleast, a fifth order accurate finite difference technique 
even if the spectral convergence is not considered. This results in a total number of 211 
points in the streamwise direction. A second-order Adams-Bashforth finite difference 
scheme was utilized for temporal discretization. With the computationally affordable 
211 x 65 grid resolution available, accurate simulations with Re = Pe = 200 (where these 
normalized quantities are defined based on the velocity difference across the layer, and 
the initial vorticity thickness) are possible [8]. These values are somewhat smaller than 
those of a fully developed laboratory turbulent shear layer, but are sufficiently high to 
promote the growth of instability waves. 

In order to visualize the flow evolution, a time series of DNS generated snapshots of the 
instantaneous normalized product concentration are shown in Fig. 1. This figure depicts 
the way in which the small amplitude perturbations manifest themselves in the formation 
of large scale coherent vortices downstream of the splitter plate. The forcing associated 
with the most unstable mode alone would cause the initial roll-up of these vortices. The 
perturbations corresponding to the first subharmonic mode result in a second roll-up in 
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the form of merging neighboring vortices, and the presence of the second subharmonic 
generates a third roll-up (second merging) at a region near the outflow. The roll-up and 
the subsequent pairings of these vortical structures engulf the unreacted species from 
the two streams into the mixing zone. The dynamics of the vortical evolution in this 
process reflect the two main mechanisms of mixing enhancement. Large concentrations of 
vorticity, that are approximately of one sign, bring the fluid from the two free streams into 
the large scale structures. Within these structures, the fluid elements are subject to further 
straining, and diffusion between the two fluids results in final mixing at the molecular 
scale. With these processes, along with the assumption of infinitely fast chemistry, the 
rate of product formation is naturally enhanced by the dynamics of vortical evolution. 

To demonstrate the capability of the model based on the Beta density, the cross stream 
variation of the product concentration and those of the negative of the unmixedness 
are presented at two streamwise locations in Figs. 2-3. Also the streamwise variation 
of the total product, T v , is presented in Fig. 4. In all these figures, the comparisons 
between the model predictions and the DNS results are noteworthy. The agreements 
demonstrated by these comparisons provide a reasonable justification for recommending 
Beta density as working relations in routine engineering predictions. However, these 
relations are advocated only for statisticanl predictions of low order moments. The 
physics of the mixing phenomenon in a spatially developing shear flow is far too complex 
to be completely described by the Beta density. The highly intermittent nature of the 
mixing process exhibits features that cannot be fully captured by this density [8]. To 
demonstrate this point quantitatively, in Figs. 5-6 the cross stream variations of the third 
and fourth centralized moments of the Shvab Zeldovich variable calculated from the /? 
density are compared with those generated by DNS. These figures indicate that while the 
trend is somewhat similar in the two cases, the extent of agreement is not as good as those 
of lower level statistical quantities (Figs. 2-4). This is primarily due to the mechanism of 
large scale entrainment in the shear flow, in that there are always unmixed fluids present 
within the core of the vortical structures. This results in a slight trimodal behavior in the 
probability distributions generated by DNS and observed experimentally [4, 9, 10, 11], 
This trimodal behavior cannot be captured by the /? density which is at best suitable for 
producing bimodal distributions. 
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4 Appendix II 


This appendix has been provided by Mr. Craig Steinberger and Mr. Thomas Vidoni. Mr. 
Steinberger is currently being partially supported under this grant. 


Problem Description 


The flow under investigation consists of a spatially developing, two-dimensional, planar 
jet flame. The fuel, F, discharges from the inner higher speed jet to a co-flowing lower 
velocity oxidizer, 0, stream. The flow evolves spatially in the streamwise direction 

(x) , and impermeable boundary conditions are imposed in the cross-stream direction 

(y) . The computational procedure involved in DNS is the same as that employed in 
[12]. This involves the discretization of the complete Navier-Stokes equations as well 
as the applicable Fickian chemical species conservation equations. This discretization 
is based on a two- four (second order accurate in time, and fourth-order truncation in 
space) compact parameter finite difference scheme. This scheme has been developed 
and made fully operational in [13] and is utilized in the context of the SPARK computer 
code developed in [14]. A generalized single step Arrhenius exothermic chemistry model 
of the form F + O — * Products + Heat is assumed. This kinetic mechanism is simple 
enough to allow an efficient treatment via DNS, and yet it is capable of capturing some 

" of the important non-equilibrium effects associated with combustion. All the species 
involved in this reaction are assumed to be thermodynamically identical, and the values 
of the physical transport parameters are assumed to be invariant with temperature. The 
flow field is initialized with a top hat streamwise velocity profile with a low amplitude 
random forcing at the jet entrance. The magnitudes of the forcing frequency and the 
phase shifts between the different components of the forcing eigenfunctions are selected 
randomly from a seed with a top hat probability frequency and with specified values of 
the mean and the rms. This random specification of the frequency and the phase shift is 
intended to mimic the influence of turbulence in an otherwise deterministic simulation. 
The amplitudes of the perturbation are chosen so that the maximum value of the turbulent 
intensity at the inflow is no more than 5%. 
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The parameters influencing the rate of reactant conversion in this setting are the 
Damkohler number, the Zeldovich number, the Reynolds number, the heat release 
parameter, the convective Mach number and the Prandtl/Schmidt numbers. Based on 
our earlier findings [15, 12], in the range of parameters considered the Reynolds number 
does not have a significant effect. Therefore, the influences of this parameter and those 
of the Peclet and the Schmidt numbers are not investigated. The magnitudes of the 
Zeldovich number and the heat release parameter are held constant in all the simulations. 
In the computational formulation, the full compressible form of the transport equations 
are considered. However, the range of the convective Mach numbers considered is not 
very large; therefore, many of the physical issues associated with high speed combustion 

[14. 16. 12] are not addressed. 

With available computational resources, numerical experiments with a resolution of 
245 (in x) xl65 (in y) mesh are possible. This resolution is sufficient provided that 
an appropriate grid generation procedure with a non-uniform mesh and with a heavy 
concentration of grid points near the regions of maximum instantaneous shear is utilized 

[14. 15. 12] . With this configuration, it is possible to perform simulations with moderate 
values of the Reynolds, Peclet, Damkohler and Zeldovich numbers within a domain 
13jD > x > 0, 3|D > y > -3|D, where D denotes the width of the jet at the inflow. 
With the numerical resolution attained, it was possible to perform simulations with a 
Reynolds number of 3 x 10 4 based on the velocity difference between the two jet streams 
at the inflow ( Uf - Uo ), the jet width (D), and the density at the inlet. The magnitudes of 
the Prandtl and the Schmidt numbers were assumed to be equal to unity, resulting in an 
identical value for the Reynolds and the Peclet numbers. The velocity ratio between the 
two streams is ^ A complete typical simulation requires about 6 hours of CPU time 
on a CRAY-2 supercomputer. 

In order to show the effects of exothermicity, simulations are performed of both a heat 
releasing (flame D) and a non-heat releasing flame (Flame A). A comparison between 
the product thickness of this two flames is made in Fig. 6. This figure indicates a 
higher rate of product formation for the case of exothermic reaction. This observation 
is not consistent with any of the previous DNS [17, 12] or experimental [18, 19] results. 
Based on our observation, we recommend that the effects of exothermicity be assessed 
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by means of laboratory measurements. These measurements must involve a reacting 
system whereby the rate of reaction conversion is temperature dependent and in which 
the large scale mixing intensity is not significantly affected by the heat release. It is 
important to note that the issue would not be resolved by a mere comparison of the 
mixing characteristics between a reacting and a non-reacting layer. An appropriate 
analysis requires the consideration of two reacting systems with the same hydrodynamic 
characteristics but with different magnitudes of the enthalpy of combustion. The extent 
of validity of our conclusion under more complex chemistry models can be determined 
by these experiments. 
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Figure Captions 


Figure 1. Time sequence of contour plots of instantaneous product concentration. 

Figure 2. Cross stream variation of mean product concentration (< C p >) from the DNS 
and the model; (a) x/8 = 11, (b) x/8 = 22. 

Figure 3. Cross stream variation of the negative of the unmixedness from the DNS and 
the model; (a) x/8 = 11, (b) x/8 = 22. 

Figure 4. Streamwise variation of total product ( T p ) from the DNS and the model. 

Figure 5. Cross stream variations of the third moment (7/3) of the Shvab-Zeldovich 
variable from the DNS and the model; (a) x/8 — 11, (b) x/8 = 22. 

Figure 6. Cross stream variations of fourth moment (^4) of the Shvab-Zeldovich variable 
from the DNS and the model; (a) x/8 = 11, (b) x/8 = 22. 

Figure 7. The streamwise variation of the product thickness under both non-heat releasing 
(Flame A) and heat releasing (Flame D) conditions. 
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